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Abstract 

Three dimensional stress analysis of spiral bevel gears in 
mesh using the finite element method is presented. A 
finite element model is generated by solving equations that 
identify tooth surface coordinates. Contact is simulated 
by the automatic generation of non-penetration constraints. 
This method is compared to a finite element contact 
analysis conducted with gap elements. 

Introduction 

Spiral bevel gears are used to transmit power between 
intersecting shafts. In helicopter transmissions these gears 
operate at very high speeds and transmit substantial 
power. To reduce vibration, noise and improve life, the 
stresses and deflections during meshing must be 
understood. 

The finite element method (FEM) is particularly well 
suited to study tooth deflections, contact stresses and tooth 
fillet stresses. Much effort has focused on predicting 
stresses in gears with FEM. Most of this research, to 
date, has been done on parallel axis gears with two 
dimensional geometry. Because of the difficulty of 
correctly identifying the three dimensional geometry, only 
a few researchers have investigated spiral bevel gears with 
FEM [1,2]. 

Previously a single spiral bevel tooth has been analyzed 
with FEM [3], however in that analysis a contact stress 
distribution was assumed. To model the correct contact 
load distribution requires a gear and pinion loaded in 
mesh. In reference [4], a gearset in mesh was analyzed 
using gap elements to model the contact zones. 


The research reported in this study will utilize the 
numerical solution for spiral bevel surface geometry with 
the automatic generation of non-penetration constraints to 
model the contacts during meshing. Pinion tooth and gear 
tooth surfaces will be developed based on the 
manufacturing kinematics. A multi-tooth model (4 gear 
and 3 pinion teeth) is created by replicating pinion and 
gear teeth in space with a solid modeler [10]. The results 
for the gap element model are compared to the non- 
penetration constraint model. 

Tooth Surface Coordinate Equations 

The system of equations, required to define the coordi- 
nates of spiral bevel gear tooth surfaces, were derived in 
reference 3, and are briefly summarized here. 

The first equation (equation of meshing) is based on the 
kinematics of manufacture and the machine tool settings. 
The equation of meshing requires the relative velocity 
between a point on the cutting surface and the same point 
on the pinion being cut must be perpendicular to the 
cutting surface normal. 

n • V = 0 

where, 

n is the normal vector from the cutter surface and 

V is the relative velocity between the cutter and the 
workpiece surfaces at the specified location. This 
equation is developed in terms of machine tool parameters 
U, 0 and where, U and 0 (length and rotational 
orientations) locate a point on the cutting head cone and 
4> c is the rotated orientation of the cutter as it swings on 
the cradle. 
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(3) 


Since the kinematic motion of cutting a gear is equivalent 
to the cutting head meshing with a simulated crown gear, 
an equation of meshing can be written in terms of a point 
on the cutting head i.e, in terms of U, 6 and ^ c . The 
equation of meshing for straight sided cutters with a 
constant ratio of roll between the cutter and the workpiece 
is given by (ref, 5): 

(U -r cotijr cos ^r) cosy sint 
+ S(w w - siny) costy sin0 
t cosy siny sin(q - (1) 

± E m (cosy siny + siny cosy cost) 

- L m siny cos^ sint = 0 

The upper and lower signs are for left and right hand 
gears respectively* The following machine tool settings 
are defined: 

\p cutting blade angle 
r (0Tq 
q cradle angle 
y root angle of workpiece 
E* machining offset 

L* vector sum of change of machine center to back 
and the sliding base 

Toot the relationship between the cradle and the 

workpiece for a constant ratio of roll 
U generating cone surface coordinate 
S radial location of cutting head in coordinate 
system S a 

r radius of generating cone surface 
Equation (1) is equivalent to: 

A (u, e, 40 = o ( 2 ) 

Because there are three unknowns U, 0, and </> c ; three 
equations must be developed to solve for the surface 
coordinates of a spiral bevel gear. The three parameters 
U, 0 and 4 C a*© defined relative to the cutting head and 
cradle coordinate systems (S c and SJ respectively. These 
parameters can be transformed through a series of 
coordinate transformations to a coordinate system attached 
to the workpiece [5]. Or U, 0 and <f> e can be mapped 
into X*, Y w , Z* in coordinate system ^ attached to the 
workpiece. These transformations, used in conjunction 
with two other geometric requirements, give the two 
additional equations. 

The correct U, 0 and <t> c that solves the equation of 
meshing, must also, upon transformation to the workpiece 
coordinate system S*, result in a axial coordinate Z* that 
matches with the value of Z found by projecting the tooth 


face in the XY plane. 

Z„ - Z = 0 

This equation along with the correct coordinate trans- 
formations (see Equation 10) result in a second equation 
of the form: 

f 2 (u, e, ) = 0 (4) 

A similar requirement for the radial location of a point on 
the workpiece results in the following: 

r - J&l * Y$ = 0 (5) 

These projections are as shown in figure 1. The 
appropriate coordinate transformations (see Equation 10) 
will convert equation 6 into a function of U, 0 and <j> c 

&u t e, a » o (6) 

Equations (2), (4), (6) form a system of nonlinear 
equations necessary to define a point on the tooth surface. 

Solution Technique 

Equations (2), (4), (6) form a nonlinear system of equat- 
ions that do not have a closed form solution. They are 
solved using Newton’s method (ref. 8). 

An initial guess U # , 0«, is used to the start iterative 
solution procedure. Newton’s method is used to 
determine subsequent values of the updated vector (U k , 

0k* ^ck)' 


U t 




y 1 , 
r *-l 

0* 

= 

e*-. 

+ 

yL 

«>* 
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where the vector Y is the solution of: 

The 3x3 matrix in the preceding equation is the Jacobian 
matrix and must be inverted each iteration to solve for the 


dfiiU*' 1 ) ^(e* -1 ) ^r i (4»e' 1 ) 

dU 06 04 » e 
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Y vector. The equation of meshing, function f v is 
numerically differentiated directly to find the terms for the 
Jacobian matrix. Function f 2 and f 3 cannot be directly 
differentiated with respect to U, 0 and 0 C . After each 
iteration U 1 * 4 , 0 M , (in the cutting head coordinate 
system) are transformed into the workpiece coordinate 
system, S w , with the series of coordinate transformations 
as given in Equations 9 and 10. (The right side of 
equation 9 locates a point on the cutting head cone as 
shown in figure 2.) 


x w 


r cotijr - U cos\|r 

Y w 

* [MJ 

U sini|r sin0 

2 W 

U smijr cos0 

1 


1 


where, 

[*_]-[ (10) 

The matrix [M*J above represents a series of 
homogeneous coordinate transformations from a 
coordinate on the cutting head onto a coordinate system 
on the workpiece [5]. 

Functions f 2 and f 3 are evaluated by starting with an initial 
U k , 0 k and performing the transformations in 
Equation (10) and evaluating Equations (4) and (6). The 
numerical differentiation of f 2 and f 3 is performed by 
transforming U k +inc, 0 k +inc, ^ k +inc (where inc is a 
small increment appropriate for numerical differentiation) 
into X^+inc, Y w +inc, Z^inc. Equations (4) and (6) 
are then used to evaluate the numerical differentiation. 
Function f v f 2 , f 3 and their partial derivatives are required 
for the Jacobian matrix and are updated each iteration. 
The iteration procedure continues until the Y vector is less 
then a predetermined tolerance. This completes the 
solution technique for a single point on the spiral bevel 
gear tooth surface. In this way the coordinates of the 
surface of the tooth are defined. This solution technique 
is repeated for each of the four surfaces; gear convex, 
gear concave, pinion convex and pinion concave. 

Since all four surfaces are generated independently, 
additional matrix transformations are required to obtain 
die correct tooth thickness. The concave surfaces are 
fixed on each tooth. The convex surfaces are rotated as 
required. The angle of rotation is obtained by matching 
the tooth top land thickness with the desired value. 


Gear and Pinion Orientations Required for Meshing 

When generating the pinion and gear surface as described 
above, the pinion cone and gear cone apex will meet at 
the same point. This point is the origin of the fixed 
coordinate system attached to the workpiece being 
generated. To place the gear and pinion in mesh with 
each other, rotations described in the following example 
are required (ref. 6). 

1. The gear tooth surface points are rotated by 
360/N t + 180 CW about the global Z*, axis, for this 
example, the rotation is 190 degrees. 

2. The pinion is rotated by 90 deg CCW about the 
global Y axis. 

3. Because the four surfaces are defined 
independently, their orientation is random with 
respect to meshing. The physical location of the 
gear and pinion after rotations described above 
correspond to the gear and pinion overlapping. 
To correct this condition the pinion is rotated CW 
about its axis of rotation Z ^ until surface contact 
occurs. For the example used in this study, the 
rotation was 3.56 deg. 

To make a complete pinion, the generated pinion tooth 
was copied and rotated 12 times and the generated gear 
tooth was copied and rotated 36 times with the aid of a 
solid modelling program as shown in figure 3. 

Finite Element Analysis of Spiral Bevel Gears 

Recent finite element analysis research on spiral bevel 
gears has attempted to solve the contact stress distribution 
in a multi-tooth model (4 gear and 3 pinion teeth) (ref. 4). 
The tooth pair contact zones in reference 3 were modeled 
with gap elements. In the current study, use of software 
with automated contact options is used. This is done to 
avoid certain limitations in the use of gap elements, such 
as: 

(i) Difficulties in connecting the gap elements with 
proper orientations. 

(ii) The orientation of the contact plane remains 
unchanged during deflection. 

(iii) Difficult to accurately select the properties such 
as appropriate open/closed stiffness values, 
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selection of the stiffness matrix update strategy 
and efficient problem restarts. 

New advances in the state of the art for FEM provide 
deformable body against deformable body penetration 
algorithms which can be used to establish the nonlinear 
boundary conditions for contact problems. MARC (ref. 
9) is one such FEM package software which uses this 
algorithm to automatically detect contact. 

Automated Contact Analysis 

Situations in which contact occurs are common to many 
different nonlinear applications. Contact forms a 
distinctive and important subset to the category of 
changing-status nonlinearities. Contact, by its very nature, 
is a nonlinear problem. During contact, both the forces 
transmitted across the area and the area of contact change. 
The force-displacement relationship thus become 
nonlinear. Usually, the transmitted load is in the normal 
direction. 

The program used MARC (ref. 9), is a general purpose 
computer program for linear and nonlinear stress and heat 
transfer analysis. This program is capable of solving 
problems with nonlinearities that occur due to material 
properties, large deformations, or boundary conditions. 
In general, the solution of nonlinear FEM problems 
requires incremental solution schemes and sometimes 
requires iterations within each load/time increment to 
satisfy equilibrium conditions at the end of each step. 

The FEM program used has a folly automatic contact 
option which enables the analysis between finite element 
bodies without the use of any special gap or contact 
elements. The contact option was originally designed for 
analysis of manufacturing processes such as forging or 
sheet metal forming, but its capabilities have been 
expanded to meet other analysis requirements. The work 
presented here utilized the program of reference 9 running 
on a supercomputer. 

Contact between the bodies is handled by imposing non- 
penetration constraints. The non-penetration constraint, 
as shown in figure 4 and figure 5 is: 

U A • n < D 

where: U A is the displacement of node A, 
n is normal to the contacted surface, 

D is the contact tolerance 


Solver constraints are used to impose the non-penetration 
constraint, and a very efficient surface contact algorithm 
which allows the user to simulate general 2D/3D multi- 
body contact. Both "deformable-to-rigid" and 
"deformable-to-deformable" contact situations are 
allowed. The user needs to only identify which bodies 
might come in contact during the analysis. The bodies 
can be either rigid or deformable, and the algorithm 
tracks variable contact conditions automatically. Thus, 
the user no longer needs to worry about the location and 
open/close status checks of "gap elements". Automatic, 
in this context, means that user interaction is not required 
in treating multibody contact and friction, and the 
program has automated the imposition of non-penetration 
constraints. 

Real world contact problems between rigid and/or 
deformable bodies are actually 3D in nature. To solve 
such contact problems, one needs to define bodies and 
their boundary surfaces. 

Deformable bodies are defined by the elements of which 
they are made. Once all the boundary nodes for a 
deformable body are determined, 4 point patches are 
automatically created, which are constantly updated with 
the body deformation. Contact is determined between a 
node and all body profiles, deformable or rigid. 

The user must define bodies so that their boundary sur- 
faces can be established. Deformable bodies are defined 
by a list of finite elements. Because the contact boundary 
conditions are a function of the applied load, the analysis 
must be carried out incrementally. Within each load or 
time increment of an analysis, additional iterations may be 
required to stabilize the contact zone. Contact problems 
involve two important aspects: 

(i) the opening and closing of the gap between bodies 

(ii) friction between the contacting surfaces. 

The MARC program establishes a hierarchy between the 
bodies so that at a given contact interface, one body is the 
contactor and the other body is the contacted. The set of 
nodes on the boundary surface of a contactor are 
candidate nodes for contact. The boundary surface of a 
contacted body is defined by a set of geometrical entities. 
A user specified contact tolerance is used to determine the 
body separation distance which determines whether two 
bodies are considered to be in contact with each other. 
The contactor’s boundary nodes are prevented from 
penetrating the surface of the contacted body by imposing 
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solver constraints. For contact between two deformable 
bodies, MARC applies multipoint constraints in the form 
of ties. The ties link the motion of one node on the 
contactor body to two adjacent nodes on the contacted 
body. During each iteration as nodes enter and leave 
contact or slide between adjacent node segments on 
contacted bodies, a bandwidth optimization is performed 
to reduce the computer processing time required. 

A static analysis of two bodies that are not initially 
connected poses special problems with a FEM, if one of 
the bodies has a rigid body motion component If, at any 
time, the two bodies are disconnected then the stiffness 
matrix would become singular and unsolvable (in a static 
analyses). In order to overcome this difficulty, the two 
gears are connected with weak springs. The spring 
stiffness must be negligible compared to the material 
stiffness. 

Model descriptions 

Pinion and gear design data and generation machine 
settings are as shown in Tables I and II. This information 
is used as input into a computer program that solves the 
3 nonlinear equations (previously described) used to 
define the points on the pinion and gear surface. The 
program generates a NxM mesh on the tooth profile of a 
spiral bevel gearset [7]. After one gear and one pinion 
tooth was made, the geometric modeler [10] was used to 
rotate the teeth to create 4 gear teeth and 3 pinion teeth in 
mesh. Eight noded, isoparametric 3D brick elements 
were used. The seven tooth model is as shown in figure 
6. The seven tooth model consisted of 8793 elements and 
11261 nodes. 

Loading and Boundary Conditions 

The boundary conditions for the seven tooth model are 
also shown in figure 6. The boundary conditions are 
applied such that the gear teeth are made to pivot about a 
fixed point along the axis of rotation for the gear on the 
Z axis. The nodes where the forces are applied are in the 
hub of the gear. Fixed displacement boundary conditions 
are applied at 8 nodal points in the Z direction only for 
the gear and in all directions for the pinion as shown in 
figure 6. Since this is a static problem involving two 
bodies (the pinion and the gear) in contact, as described 
in a previous chapter, weak springs are added to prevent 
the rotational rigid body modes for the gear. Eight 
springs are added away from the contact region. The 
springs connect the comer nodes of the pinion with the 
comer nodes of the gear on the faces where contact 


occurs. The stiffness of the springs are 100 lbs/in. 
This is insignificant when compared to the tooth contact 
stiffness and therefore does not effect the overall solution. 
The maximum torque for the gear mesh studied was 9508 
in. •lbs. on the gear. This torque is applied as a 
concentrated force with a moment arm on the gear hub. 
This concentrated force for the seven tooth model was 
4714 lbs. with a moment arm of 2.017 inches. The 
force is applied incrementally for convergence to occur. 

Gap Element and Automated 

Contact Analysis Comparison 

Contact stresses on spiral bevel gears were studied by 
researchers with gap elements in references 11 and 12 on 
a similar seven tooth model, comparison of the results 
with automated contact analysis is presented. Both 
models contained the same mesh density, boundary 
conditions, material properties and loading. The nodal 
stress results of pinion tooth #1 obtained from automated 
contact analysis and gap elements are as shown in figures 
7 and 8, respectively. In both analyses the highest 
concentrated stress value occur at the same node. 
Comparison of the results of these two runs are as 
follows. 



Gap Element 

Automated 
Contact Analysis 

Min nodal principal 
stress 

-296,410 (psi) 

-291,503 (psi) 

CPU time (approx) 

30 min. 

80 min. 

Elemental principal 
stress 

-84,761 (psi) 

-113,577 (p*i) 

Gap element closed or 
nodal points with contact 

4 

8 

No. of iterations 

4 

8 


The nodal stresses were in very good agreement for both 
the gap element FEM and the automated contact FEM 
model. Elemental and nodal stresses differed greatly for 
both models indicating a need for mesh refinement. The 
gap element model, as expected, used less cpu time. 
More information is given by virtue of locating and 
defining the gap elements. This reduced the 
computational time required to search for and calculate 
the contact area and stresses. However, as noted earlier, 
it is much more difficult to calculate the proper 
orientation and location of the gap elements. It is more 
convenient to analyze contact with the automatic contact 
FEM program. 
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The number of contact nodes at the contact region in the 
■ntnmatffd contact FEM program was higher than 
identified by die gap element FEM program. With the 
pinion considered body 1 (contactor body) and the gear 
considered body 2 (contacted body), eight nodes contacted 
as shown above. With body 1 and body 2 switched, 16 
nodes were found to have contact. Presumably this sort 
of discrepancy occurred because the mesh was too coarse. 

CONCLUSIONS 

Three-dimensional contact analysis of spiral bevel gears 
in mesh was performed. Four gear teeth and three pinion 
teeth, generated by solving equations based on gear 
manufacturing kinematics, were used to model contact. 
An automatic contact analysis algorithm is utilized with 
the finite element method. The contact algorithm searches 
for and finds contact without using gap elements. 

Comparison of nodal and elemental principal stresses is 
favorable with a previous FEM that utilized gap elements. 
However, for both models the nodal and elemental 
stresses were significantly different. Also, both models 
had large stress gradients at the contact point. This 
indic ates that a more refined finite element mesh is 
required. 

The contact algorithm model takes more cpu time 
compared to the gap element model. The contact 
algorithm has to find where contact occurs, with gap 
dements contact is defined to occur at the gap elements. 

Although the gap element modal ran with less cpu time, 
use of the contact algorithm is more convenient. The 
contact algorithm is easier because gap elements do not 
have to be mathematically defined in space with proper 
orientation. 
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TABLE I: PINION AND GEAR DESIGN DATA 


| PINION 

GEAR 

| — — 

Number of teeth pinion 

12 

36 

Dedundum angle, degree 

1.5666 

3.8833 

Addendum angle, degree 

3.8833 

1.5666 

Pitch angle, degree 

18.4333 

71.566 

Shaft angle, degree 

90.0 

90.0 

Mean spiral angle, degree 

35.0 

35.0 

Face width, mm (in) 

25.4 (1.0) 

25.4 (1.0) 

Mean cone distance, mm (in) 

81.05 (3.19) 

81.05 (3.19) 

Inside radius of blank, mm (in) 

5.3 (0.6094) 

3.0 (.3449) 

Top land thickness, mm (in) 

2.032 (0.080) 

2.489 (0.098) 

| Clearance, mm (in) 

0.762 (0.030) 

0.929 (0.0366) 


TABLE H: GENERATION MACHINE SETTINGS 
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Figure 5. Contact depiction 
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Figure 6. Seven tooth model with boundary conditions 
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Figure 7. Nodal stress result on pinion obtained from 

automated contact analysis in seven tooth model 
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Figure 8. Nodal Stress Result on Pinion Obtained From the Gap 
Element Solution and Seven Tooth Model 
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